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Abstract: We propose an estimation approach to analyse correlated functional 

data which are observed on unequal grids or even sparsely. The model we use is a 
functional linear mixed model, a functional analogue of the linear mixed model. Esti¬ 
mation is based on dimension reduction via functional principal component analysis 
and on mixed model methodology. Our procedure allows the decomposition of the 
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variability in the data as well as the estimation of mean effects of interest and bor¬ 
rows strength across curves. Confidence bands for mean effects can be constructed 
conditional on estimated principal components. We provide R-code implementing our 
approach. The method is motivated by and applied to data from speech production 
research. 


Key words: Functional additive models, functional data, functional principal com¬ 
ponent analysis, penalized splines, speech production. 


1 Introduction 

Advancements in technology allow today’s scientists to collect an increasing amount 
of data consisting of functional observations rather than single data points. Most 
methods in functional data analysis (fda) (see e.g. Ramsay and Silverman, 2005) as¬ 
sume that observations are a) independent and/or b) observed at a typically large 
number of the same (equidistant) observation points across curves. 

Linguistic research is only one of numerous fields in which the data often do not 
meet these strong requirements. Our motivating data come from a speech production 
study (Pouplier et ah, 2014) on assimilation, the phenomenon that the articulation of 
consonants becomes more alike when they appear next to each other or across words. 
The data consist of the acoustic signals of nine speakers reading the same 16 target 
words including the consonants of interest each five times. Due to the repeated mea¬ 
surements for speakers and for target words, the data have a crossed design structure. 
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Moreover, number and location of the observation points differ between the observed 
curves. 

We propose a model and an estimation approach that extend existing methods by ac¬ 
counting for both a) correlation between functional data and for b) irregular spacing 
of - possibly very few - observation points. The model is a functional analogue of 
the linear mixed model (LMM), which is frequently used to analyse scalar correlated 
data. 

For dimension reduction, we use functional principal component analysis (FPCA, see 
e.g. Ramsay and Silverman, 2005) to extract the dominant modes of variation in the 
data. The functional random effects are expanded in bases of eigenfunctions of their 
auto-covariances, which we estimate beforehand using a method of moments estima¬ 
tor. FPCA has become a key tool in fda as it yields a parsimonious representation of 
the data. It is attractive as the eigenfunction bases are estimated from the data and 
have optimal approximation properties for a fixed number of basis functions. It also 
allows for an explicit decomposition of the variability in the data. 

We propose two ways of predicting the eigenfunction weights. We either compute 
them directly as empirical best linear unbiased predictors (EBLUPs) of the resulting 
LMM or we alternatively embed our model in the more general framework of func¬ 
tional additive mixed models (FAMMs) introduced by Scheipl et al. (2015). The first 
approach is straight forward and computationally more efficient; it does not require 
additional estimation steps as a plug-in estimate is used. The latter has the advantage 
that all model components are estimated/predicted in one framework and it allows 
for approximate statistical inference conditional on the FPCA. 

Previous work for dependent functional data differ in their generality and restrictions 
on the sampling grid, but only one approach can cope with correlated functional data 
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observed on unequal grids and with complex correlation structures. Brumback and 
Rice (1998) consider a smoothing spline-based method for nested or crossed curves, 
which are modelled as fixed effect curves. They allow for missing observations in 
equal grids but do not consider any covariate effects. A Bayesian wavelet-based func¬ 
tional mixed model approach is introduced by Morris et al. (2003) and extended by 
Morris and Carroll (2006), Morris et al. (2006), and subsequent work by this group. 
While this approach is quite general in the possible functional random effects struc¬ 
ture, it assumes regular and equal grids with at most a small proportion of missings 
and a reasonable number of completely observed curves. Di et al. (2009), Greven 
et al. (2010), and Shou et al. (2015) consider functional linear mixed models with 
a functional random intercept, with a functional random intercept and slope, and 
with nested and crossed functional random intercepts, respectively. While following 
a similar approach to estimation for these models without the options for approxi¬ 
mate inference we provide, all three are restricted to data sampled on a fine grid. Di 
et al. (2014) extend the random intercept model of Di et al. (2009) to sparse func¬ 
tional data; the correlation structure, however, remains less general than ours. Also 
motivated by an application from linguistics, Aston et al. (2010) perform an FPCA 
on all curves ignoring the correlation structure and then use the functional principal 
component (FPC) weights as the response variables in an LMM with random effects 
for speakers and words. Only linear effects of scalar covariates are considered, FPC 
bases are restricted to be the same for all latent processes, and it is assumed that the 
data are sampled on a common grid. Brockhaus et al. (2015) propose a unified class 
for functional regression models including group-specific functional effects, which are 
represented as linear array models and estimated using boosting. The array structure 
requires common grids. Other approaches concentrate specifically on spatially corre- 
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lated functional data on equal grids, as e.g. Staicu et al. (2010). A very flexible class 
of functional response models has been developed by Scheipl et al. (2015), allowing 
for multiple partially nested and crossed functional random effects with flexible cor¬ 
relation structures. Both spline-based and FPC-based representations incorporated 
in a mixed model framework are considered. They allow for densely and sparsely 
sampled data, but in the case of the FPC-based representation, they assume that ap¬ 
propriate FPC estimates are available. We take advantage of this general framework 
and use our estimated FPCs to obtain improved estimates and approximate point- 
wise confidence bands for the mean and covariate effects. In addition to providing an 
interpretable variance decomposition, our FPC-based approach reduces computation 
time by orders of magnitude compared to the spline-based estimates from Scheipl 
et al. (2015) (compare Section 5), allowing the analysis of realistically sized data in 
practice. 

A number of approaches allow for irregularly or sparsely sampled functional data 
but assume that curves are independent. Guo (2002, 2004) first introduce the term 
functional mixed effects models for their model. The model does not capture between- 
function correlation as only curve-level random effect functions are included which 
are modelled using smoothing splines. The approach is not restricted to regularly 
sampled grid data. Chen and Wang (2011) propose a spline-based approach that is 
suitable for sparsely sampled data, but similar to Guo (2002, 2004) they only consider 
curve-level random effects. James et al. (2000), Yao et al. (2005), and Peng and Paul 
(2009) among others propose FPCA approaches for sparsely observed functional data 
with uncorrelated curves. 

For an extensive overview and further references for functional regression approaches, 
including functional response regression, see Morris (2015). 
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The remainder of the paper is organized as follows. Section 2 introduces the gen¬ 
eral functional linear mixed model and presents an important special case which is 
used to analyse the motivating assimilation data. Section 3 develops our estimation 
framework. Our method is evaluated in an application to the assimilation data and 
in simulations in Sections 4 and 5, respectively. Section 6 closes with a discussion 
and outlook. Theoretical results and supplementary material including estimation 
details as well as additional results for application and simulations are available in 
the appendix. 


2 Functional linear mixed models 

2.1 The general model 

The general functional linear mixed model (FLMM) is given by 

Yi(t) = n(t,Xi) + zJU(t) + Ei(t) + £i(t), i = 1,... ,n, (2.1) 

where Y l (t) is the square-integrable functional response observed at arguments t in 
T, a closed interval in M, and n is the number of curves. fi{t, ®j) is a fixed main effect 
surface dependent on a vector of known covariates x t . To account for the functional 
nature of the Y)(f), the random effects of an LMM are replaced by a vector-valued ran¬ 
dom process U(t). Zi is a known covariate vector. E t {t) is a curve-specific deviation 
in form of a smooth residual curve. We assume that there is white noise measurement 
error denoted by c*(f) with variance a 2 that captures random uncorrelated variation 
within each curve. Note that if needed, the error variance may also vary across t, 
cr 2 (f). We further assume that U(t), £)(£), and £i(t), i = 1 , ...,n, are zero-mean, 
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square-integrable, mutually uncorrelated random processes. Therefore, each of the 
q components of U(t) has an auto-covariance operator K Uj (t,t'), j = 1 ,q, and 
cross-covariance operators K u ^ k (t, t'), j,k = 1,,q, some of which might be zero 
for uncorrelated functional random effects. Ei(t) has an auto-covariance operator 
K E (t,t') = Cov [Ei(t), Ei(t')]. In the following, mean, auto-covariances, and thus 
also the eigenfunctions are assumed to be smooth in t. 

p(t,Xi) is an additive function of t and x*. For example, it can be constant in t, 
fi(t,Xi) = p(xi), or additive in t and Xj, p(t,Xi) = p\(t) + /i 2 (x,). Another spe¬ 
cial case is when all Xu, ..., Xi P in x, act as index-varying coefficients, pit, x,) = 
fo(t) + fi(t)xn + ... + fp(t)xi p , with unknown smooth functions fo(-), ■ ■., f p (-)- 

2.2 Special case: the FLMM for a crossed design 

For our application in speech production research (Section 4), we use an FLMM with 
a crossed design structure to account for correlation between measurements of the 
same speaker and between measurements of the same target word. 


Yijh{t) — p(t, x ijh ) + Bi(t ) + Cj(t) + Eijh(t ) + £ijh(t), (2-2) 

with i = 1 ,...,/ (number of speakers), j = 1,... ,J (number of target words), and 
h = 1,..., Hij (number of repetitions). Here, Yijhit ) is the hth curve for speaker i 
and target word j observed at time t. Bi(t ) and C 3 (t) are functional random inter¬ 
cepts (fRIs) for the speakers and target words, respectively. Curve-specific deviations 
are accommodated by the smooth residual term Eijh(t), which also captures interac¬ 
tions between speakers and words. We decided to not include an interaction effect 
separately based on substantive considerations and the limited sample size. £ijh(f) 
is additional white noise measurement error with variance cr 2 . We denote the auto- 
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covariance operators by K B (t,t') = Cov [Bi(t), Bi(t')], K c (t,t') = Cov [Cj(t), Cj(t')\, 
and K E (t, t') = Cov [E ijh (t), E ijh (tf)], i =, 1,..., I, j = 1,..., J, h = 1,..., H i:j . 

2.3 Irregularly and sparsely sampled functional data 

Let us now assume that for our general model (2.1) we have observed n curves on 
observation points {U i,..., Ud^ G T, i = 1,..., n. The number and the location of 
the observation points are allowed to differ from curve to curve. In the extreme, only 
one point may be observed for a curve. Moreover, the observation points of a curve 
do not have to be equally spaced. We denote realizations of the functional response 
Yi(t) at point tij by y ltlJ , j = 1,..., D t . Accordingly, we denote realizations of the 
response in (2.2) by y ijht , with t G {t ijhl ,..., t ijhD .. h }. 


3 Estimation 

We base our estimation on FPCA, which provides the dimension reduction so im¬ 
portant for functional data. In addition, FPCA allows an explicit decomposition of 
the variability in the data. Compared to other basis approaches e.g. using splines, 
the advantage of FPCA is that the eigenbases are optimal in the sense of giving the 
best approximation for a given number of basis functions. An eigen decomposition of 
the auto-covariances of U(t ) and Ei(t) based on Mercer’s theorem (Mercer, 1909) is 
used. The estimated eigenfunctions, also known as functional principal components, 
describe the main modes of variation of processes U (t) and E,(t) and the estimated 
eigenvalues quantify the amount of variability explained by the corresponding FPCs. 
To pool information across observations, which is particularly important in the case 
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of irregularly or sparsely sampled functional data, we use smoothing of the auto¬ 
covariances, cf. Yao et al. (2005) for non-correlated sparse functional data. The four 
main steps of our estimation procedure are described in the following. 

Step 1 We estimate the mean, p(t,Xi), using penalized splines based on a working 
independence assumption. 

Step 2 For each of the functional random effects, we estimate the respective auto¬ 
covariance using a smooth method of moments estimator. Subsequently, we 
evaluate each of the estimates on a pre-specified equidistant, dense grid. 

Step 3 We conduct an eigen decomposition of each of the estimated auto-covariance 
matrices. Using the Karhunen-Loeve (KL) expansion (Karhunen, 1947; Loeve, 
1945), we can represent the functional random effects using truncated bases of 
the respective estimated eigenfunctions. 

Step 4 We propose two ways of predicting the random basis weights. 

Step 1, step 3, and the first option for the prediction of the basis weights are analogous 
to the estimation proposed in Di et al. (2009), Greven et al. (2010), and Shou et al. 
(2015) for functional data sampled on an equal, fine grid and in Di et al. (2014) for a 
simpler model. For simplicity, we focus in the remainder of this section on model (2.2). 
R-code implementing our approach for model (2.2) can be provided upon request. 

3.1 Estimation of the mean function 

We estimate the mean /.i(h ayy,) based on the working independence assumption 


Yijh(t) VCijh) T ^ijh (^) 


(3.1) 


10 Jona Cederbaum et a 1. 


with i.i.d. Gaussian random variables £ijh(t). Consider the additive mean predictor 
ia(t,x ijh ) = f 0 (t) + Ep=i fp(t) x ijh p (and / J,(t,Xi ) = f 0 (t) + Ef=i f P if)x ip for model 
(2.1)). We represent the unknown, smooth functions f p (-) using splines and use a 13- 
Spline basis in the following. Let 0 denote the Kronecker product and let /u, denote 
the stacked mean of length 2D. // is then approximated by 

p 

M = E (*' ® !h) ' ^t^ P i (3.2) 

p=0 

where denotes an inflated vector of length 2D = E,: 3 /-, Api of covariate values, 
of dimension 2) x K p comprises the evaluations of the K p basis functions on the 
2D time points Ujh , and 1 kp = (1,..., 1) T of length I\ p . 6 P is a coefficient vector of 
length K p . The concrete forms of fl/j) and are given in Appendix A. We control 
the trade-off between goodness of fit and smoothness by adding a difference penalty 
(Eilers and Marx, 1996). Using the penalized splines approximation of model (3.1) 
allows us to represent the model as a scalar LMM which has the advantage that 
the smoothing parameter can be estimated as a variance component using restricted 
maximum likelihood (REML, Patterson and Thompson (1971), cf. Ruppert et al. 
(2003), sec. 4.9). 

For the practical implementation, we build on existing software and either use R- 
function (R Development Core Team, 2014) gam or function bam, both implemented 
in the R-package mgcv (Wood, 2011). The latter is especially designed for large 
datasets. Avoiding the construction of the complete design matrix leads to a much 
lower memory footprint and the possibility of parallelization also gives a considerable 
speed-up in computation time. For further details and more general mean models see 
Wood et al. (2015). In a next step, we center the data using the estimated mean, 
fi(t, Xijh)j and obtain ijijht ■ Vijht /((f, (Ey/j) ■ 
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3.2 Estimation of the auto-covariances 


The estimation of the auto-covariances is the most challenging step in the estimation 
procedure, also in terms of computational effort. We exploit the fact that for cen¬ 
tred data, the expectation of the cross products corresponds to the auto-covariance 


Cov 


Cov 


Yijh{t),Yvj'h>W) 


, which can be decomposed as follows 


= E 


Yijh(t)Yi'j'h'(t') 


(3.3) 


= + [K E (t,tj + SwSjfSui, 


with 6 xx i equal to one if x = x' and zero otherwise. (3.3) can be seen as an ad¬ 
ditive, bivariate varying coefficient model, in which the auto-covariances are the 
unknown smooth bivariate functions to be estimated, while Sm, 8jj>, 5u'5jjf5hh', and 
Su'Sjj'Shh'fitt' represent the covariates. Under a working assumption of independence 
and homoscedastic variance of the cross products, we can use each empirical prod¬ 
uct yijhtVi'j'h't ' f° r which at least i = i' or j = f to obtain smooth estimates of 
K B (t,t'), K c (t,t'), and K E {t,t') and an estimate of the error variance a 2 . The total 
number of products yijhtVi'j'h't 1 used for the estimation of the auto-covariance then 

is E',1 (Eh, E&, Diik) 2 + Ebi (E',1 E& Dijhf - E',1 Ehi (e£, D, jh )\ 

For equal Hij = H and D i:i i, = D, this would reduce to S 2 ( l /i + l /j — l /u)-, with X> 
the total number of observation points. 

For the estimation of the auto-covariances, we use bivariate tensor product P-splines 
(see e.g. Wood, 2006, sec. 4.1.8). The idea is to combine low rank marginal bases 
for each t, t' in order to obtain smooth functions of the two covariates. Then, given 
the appropriate ordering of the parameter vector, the part of the design matrix cor¬ 
responding to K x (t,t'), X e {B,C,E}, is given by the respective indicator matrix 
multiplied entry-wise by (M,*01j) ■ (1where M t x and M x denote 
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the corresponding marginal spline design matrices of rank K for covariate t and if. 
A smoothness penalty is introduced in order to avoid over-fitting. As we estimate 
auto-covariances which are symmetric by definition and therefore are naturally on 
the same scale, we choose an isotropic penalty with the penalty matrix of the form 
Sff/ = S t ® S t where S t and S t > represent the respective marginal penalty matrices 
for covariate t and t'. We use marginal B-spline bases combined with marginal differ¬ 
ence penalties. As for the mean estimation, we take advantage of the mixed model 
representation of the additive model (3.3) using REML in order to obtain estimates 
for the tensor product basis coefficients and the smoothing parameter. 

Again, we make use of existing software by applying function bam in R-package mgcv. 
Negative estimated values of a 2 are set to zero for the final estimate. Symmetry of 
the auto-covariances is ensured through the model apart from numerical inaccuracies. 
During the auto-covariance estimation, strength is borrowed across all curves. This 
can be extremely advantageous for sparse functional data when some curves only have 
very few measurements and smoothing of curves would be infeasible. 

3.3 Eigen decompositions of estimated auto-covariances 

Based on Mercer’s Theorem, the eigen decompositions of the auto-covariances are 

OO 

K x (t,f) = («'), x € {B,C,E }, 

k =1 

where > is* > ... > 0 are the respective eigenvalues, k G N. The corresponding 
eigenfunctions {(/)%,k G N}, X G {B,C,E}, form an orthonormal basis in the Hilbert 
space L 2 \T\ with respect to the scalar product 


{/. g) = J f(t)g(t)dt. 


( 3 . 4 ) 
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In practice, the smooth auto-covariances are evaluated on an equally spaced, dense 
grid of pre-specified length D. The resulting matrices are in the fol¬ 


lowing denoted as K B = K B (t d ,t d ') 


J d,d'=l,...,D 


, K c = 


K c (t d ,t d/ 


Jd,d'=l,...,D 


, and 


. We conduct an eigen decomposition of each of the 


K E = k E (t d ,t d ,; 

L J d,d'=l,...,D 

estimated auto-covariance matrices yielding estimated eigenvectors and estimated 
eigenvalues. Rescaling of the estimated eigenvectors and thus also of the estimated 
eigenvalues is necessary to ensure that the approximated eigenfunctions are orthonor¬ 
mal with respect to the scalar product (3.4). More details are given in Appendix B. 
Negative estimated eigenvalues are trimmed to guarantee positive semi-definiteness. 


Truncation of the FPCs . While in theory, there is an infinite number of eigenfunc¬ 
tions, dimension reduction and selection of the number of FPCs for each random 
process is necessary in practice. This truncation has a theoretical justification and 
can be seen as a form of penalization (see e.g. Di et al., 2009; Peng and Paul, 2009). 
Among the multiple proposals in the literature (see for an overview Greven et al., 
2010), we base our choice on the proportion of variance explained. This allows us to 
quantify the contribution of the random processes to the variation in the observed 
data. It is based on the variance decomposition of the response 

« oo oo oo 

/ Var [Y ijh {t)\ dt = J2 u k+J2 u % + J2 u k +(j2 \ r \- 
, ' r k =i k =i k =i 

The sums Ylki u k ■> v ki anf l Ylh=i u k quantify the relative importance of each 

of the three random processes. We choose principal components of decreasing impor¬ 
tance until a pre-specified level of explained variation is reached. See Appendix B for 
details and Appendix A for the derivation of the variance decomposition. 


Approximation of the functional random processes. Based on the truncation, we use 
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KL expansions to obtain parsimonious basis representations for the random processes 

n b n c n e 

B i(t) = C j(t ) = Ei ^ = ^2^3hk ( t ) k{ t )- 

k=1 k =1 /c=l 

Note that in the case of irregularly or sparsely sampled data, the observation points 
also depend on speaker i, target word j, and repetition h. Throughout this paper, 
however, we omit the additional indices for better readability. For the same reason, 
we do not emphasise that the truncation lags and eigenfunctions are estimated. By 
construction, the basis weights or FPC weights t ^ k , and £fj hk are uncorrelated 

random variables with zero mean and variance u k , k € N, X e {B, C, E}. 

For prediction of the FPC weights, we first linearly interpolate the chosen eigen¬ 
functions such that they are available on the original observation points. Due to 
the smoothness of all model components, this leads to a small error which could be 
further decreased, if desirable, by further increasing the size of the grid D. 

3.4 Prediction of the basis weights 

The basis weights give insight into the individual structure of each grouping level and 
can be used for further analyses, e.g. for classification. For a centred random process 
Xi (t), the basis weights of the KL expansion can be written as the scalar product of 
Xi (t) and (p k (t) and are often predicted using numerical integration. For irregularly 
or even sparsely sampled data, however, numerical integration would not work (well). 
Moreover, for correlated functional data contaminated with additional measurement 
error, the separation of the weights belonging to the different basis expansions re¬ 
mains unclear and ignoring the measurement error leads to biased predictions. 

These considerations motivate our two proposals for the prediction of the basis 
weights. The first is very straight-forward and computationally efficient. It general- 
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izes the conditional expectations introduced by Yao et al. (2005). The second involves 
higher computational costs but has the advantage that the mean is re-estimated in 
the same framework and that it allows for approximate statistical inference, e.g. for 
the construction of point-wise confidence bands (CBs). 

Prediction of the basis weights as EBLUPs. Using the truncated KL expansions of 
the random processes, we can approximate model (2.2) by 


n b 


N c 


N e 


Y ijh (t) = fi(t, x ijh ) + ^2 &k (t) + &k (t) + tfjhAk (t) + £ ijh(t ) (3.5) 


k =i 


k =i 


k =i 


for the discrete observation points t e {Ujhi, ■ ■ ■ , tijhD ljh }■ The resulting model (3.5) 
is a scalar LMM in which the random effects correspond to the basis weights (Di 
et ah, 2009). Once all other components are estimated, we do not need to fit the 
LMM (3.5), but can predict the basis weights as EBLUPs as derived in the following. 
Without normality assumption, they remain best linear predictors. 

Let Y denote the stacked centred response vector of length 29. Let L x e {I, J, n} and 
N x € {-/V s , N c , N e } denote the levels of the grouping variable and the truncation lag 
for process X, X £ {B,C,E}, respectively. We define £ = , with 

£ x = ^| VT ,... the stacked vector of the basis weights of length L x N x . 

Thus £ is a vector of length 91 := IN B + JN C + nN E . <I> is the joint design matrix 


of dimension 29 x 91 of the form <E = 


, where ff> s , ff> c , and dW are the 


4» s |4» c, |$ s 

respective design matrices containing the rescaled FPC estimates evaluated on the 
original observation points. The covariance matrix of £ is a diagonal matrix with 
diagonal elements corresponding to the eigenvalues of the three random processes. 
We denote its estimate in the following by G. A more detailed description of the 
form of the matrices can be found in Appendix A. Then, the EBLUP for the basis 
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weights in model (2.2) is given in the usual form by 

£ = G4> t Cov (y^ _1 Y = G$ t (& 2 1® + _1 Y (3.6) 

Equation (3.6) requires the inversion of the estimated covariance matrix of Y, which 
is of dimension 2) x 2D. For large numbers of observation points, this can be com¬ 
putationally demanding. Furthermore, when d 2 is estimated to be close to zero, the 
covariance becomes singular. Transformations with the Woodbury formula yield the 
more favourable form £ = (d 2 G + 4> T 4>j <F T Th for which the inversion is simpli¬ 
fied to that of an 91 x 91 matrix which has full rank when either d 2 is positive or 
when has full rank. In practice, when neither of these requirements is met, the 

Moore-Penrose generalized inverse is used. Note that when cr 2 is estimated to be zero, 
the EBLUP simplifies to the least-squares estimator. 

One drawback of this computationally efficient way of predicting £ is that the mean 
is estimated using a working independence assumption. This may not be statistically 
efficient and does not provide the basis for statistical inference. This motivates our 
second proposal. 

Prediction of the basis weights using FAMMs. The second option uses the fact that 
model (3.5) together with the distribution of the basis weights implied by the KL 
expansion falls into the general framework of a FAMM (Scheipl et ah, 2015) using 
suitable marginal bases and penalties. We can write model (3.5) using estimated 
eigenfunctions and -values as 

p 

Y =£ (*S ® •*’?»’’ + A )«*+£, <3.7) 

P =0 X&{B,C,E} 

with e ~ A/"(0, (J 2 I^). Y is the stacked uncentred response vector of length 2D, the 
mean is specified as in (3.2), and dW denotes an inflated 2D x L x matrix of grouping 
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indicators. The T> x N x matrix Sfr x comprises the evaluations of the N x respective 
estimated eigenfunctions on the original observation points. The concrete forms of 
the matrices x and 'F x are given in Appendix A. Adding penalties of the form 
£ xT (I^ 0 P x ) £ x , with P x = diag (pi ,..., z>^ x ) \ corresponds to the distribu¬ 
tional assumption ~ Af (0,dmg(v x ,... ,u x x )), l = 1 ,... L x , X E {B,C,E}, 
implied by the KL expansions. This set-up using linear combinations of row tensor 
product bases with an appropriate penalty falls naturally into the framework of a 
FAMM and was in fact discussed in Scheipl et al. (2015) without, however, provid¬ 
ing an approach to estimation of the eigenfunctions needed in ^f x . Mean and basis 
weights are then estimated/predicted within one framework. (3.7) is a scalar additive 
LMM, which allows to take advantage of established methods for estimation and for 
statistical inference (for more details, see Scheipl et ah, 2015). In particular, it allows 
us to construct point-wise CBs for the mean and for covariate effects. Note that the 
inference is conditional on the estimated FPCA, i.e. it does not account for the un¬ 
certainty of the estimation of the eigenfunctions and -values and for the truncation. 
In practice, we use the R-function pffr that Scheipl et al. (2015) provide in the R- 
package refundDevel (Crainiceanu et ah, 2014). Function pffr is a wrapper function 
for function gam and for related functions in the package mgcv and therefore builds 
on existing flexible and robust software. We use the pffr-constructor pcre() for FPC- 
based fRIs for the prediction of the random processes. A constraint on the functional 
random effects assures that they are centred. The resulting point-wise CBs are with a 
constraint correction (Marra and Wood, 2012). In addition to the parsimonious basis 
of eigenfunctions, this approach has the advantage of not necessitating the estima¬ 
tion of any smoothing parameters for the random processes, as the variances of the 
random weights have already been estimated and the smoothing parameter can be 
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set to one (cf. Appendix B for details). These two features lead to a drastic decrease 
in computational cost compared to spline-based prediction of the random processes, 
as is shown in our simulations in Section 5. 

The estimation quality can be further improved, if desirable, by applying the four 
estimation steps iteratively. Several possibilities are described in Appendix B, where 
also further details on the estimation and implementation can be found. 


4 Application to the speech production research data 

4.1 Background and scientific questions 

In linguistics, the term assimilation refers to the common phenomenon whereby a 
consonant becomes phonetically more like an adjacent, usually following consonant. 
For example, assimilation occurs commonly in English phrases such as “Paris show” 
in which the word-Enal /s/-sound is, in fluent speech, pronounced very similar to the 
following, word-initial /sh/-sound (Pouplier et ah, 2011). Assimilation patterns are 
conditioned by a complex interaction of perceptual, articulatory and language-specific 
factors and are therefore a central research topic in the speech sciences. In order to 
investigate assimilation in German, Pouplier et al. (2014) recorded the acoustic sig¬ 
nals of / = 9 speakers reading the same J = 16 target words, each five times. Due 
to recording errors, for some combinations only four repetitions are included in the 
data, i.e. Hij G {4,5}. The authors concentrated on variation in assimilation pat¬ 
terns for the consonants /s/, /sh/ as a function of their order (/s#sh/versus /sh#s/, 
where # denotes a word boundary), lexical stress and vowel context. Target words 
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consisted of bisyllabic noun-noun compounds. In half of the target words consonant 
/s/ is followed by word-initial /sh/, such as in the word “Callas-Schimmel”. The 
other half contains the sequence /sh#s/, e.g. “Gulasch-Symbol”. In the following, 
we will refer to the syllables containing the consonants of interest as final and initial 
target syllables (and correspondingly to final and initial target consonants). The time 
interval in which the consonants of interest appear in the utterance was cut out man¬ 
ually from the audio recording for each repetition and the time-varying signal was 
summarized in a functional index over time, varying between +1 and —1. The index 
was calculated such that for both orders it ranges from +1 for the first consonant of 
the sequence to —1 for the second consonant of the sequence (for more details, see 
Pouplier et al. (2011) and Appendix C for data pre-processing). The resulting index 
curves are displayed in Figure 1. 

A special focus lies on the asymmetry arising from the order of the consonants. We 
investigate under which conditions (order, syllable stress, vowel context) the two con¬ 
sonants assimilate and whether assimilation is symmetric with respect to the orders 
/s#sh/and /sh#s/. A common approach to the analysis of these kind of data is to 
extract curve values at pre-defined points on the time axis (e.g. 25%, 50%, 75%) which 
are subsequently used in multivariate methods (e.g. Pouplier et al., 2011). Such anal¬ 
yses fail to capture the continuous dynamic change characteristic of speech signals. 
Applying our fda-based method allows us to take into consideration the temporal dy¬ 
namics and to account for the complex correlation structure in the data which arises 
from the repeated measurements of speakers and of target words. Moreover, we can 
quantify the effect of covariates and interactions and obtain a variance decomposition. 
All utterances were recorded with the same sampling rate (32768 Hz) and then stan¬ 
dardized to a [0,1] interval as the speaking rate and hence the target consonant dura- 
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Curves for order /s#sh/ 


00 


Curves for order /sh#s/ 




04 06 

normalized time (t) 


04 06 

normalized time (t) 


GulaschSimpel 


GouacheSymbol 

GemischSalat 

GemischSalbe 


Figure 1: Index curves of the consonant assimilation data over time. Left [right]: 
Curves of order /s#sh/[/sh#s/]. Positive values approaching +1 indicate a reference 
/s/ [/sh/] acoustic pattern, while negative values approaching —1 indicate a reference 
/sh/ [/s/] acoustic pattern. 

tion differs across speakers and repetitions. After standardization, measurements are 
unequally spaced for different curves. In some data settings, landmark registration 
can be used to account for variation in time. For this application, however, regis¬ 
tration cannot replace the standardization of the time interval as different transition 
speeds between the two consonants are part of the research question of interest, i.e. of 
the assimilation process. 

4.2 A model for the consonant assimilation data 

In order to account for the repeated measurements of speakers and target words, we 
fit an FLMM with crossed fRIs, model (2.2), to the consonant assimilation data. The 
number of measurements per curve Dijh ranges from 22 to 57 with a median of 34. 
During estimation, we truncate the numbers of FPCs using a pre-specified proportion 
of explained variance of 0.95 as described in Section 3. The equidistant grid on which 
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the auto-covariances are evaluated is of length D = 100. We use cubic B-splines with 
third order difference penalties for the estimation of the mean effects and as marginal 
basis functions for the estimation of the auto-covariances. We predict the FPC weights 
using both options. As confidence bands for the covariate and interaction effects are 
of interest here, the focus lies on the second approach using the FAMM framework. 

Covariate effects . We consider four dummy-coded covariates: consonant order (or¬ 
der), stress of the final (stressl) and of the initial (stress 2 ) target syllable, which can 
be strong or weak, and vowel context (vowel), which refers to the vowels immediately 
adjacent to the target consonants and is either of the form ia or ai, e.g. Callas- 
Schimmel. Moreover, we include the interactions of the consonant order with each 
of the other three covariates. All covariates enter the mean as varying coefficients 
yielding 


A i{t, x ijh ) = f 0 (t) + fi(t) • order, + f 2 (t) • stressl, + f 3 (t) • stress2 i (4.1) 
+ / 4 (f) • vowelj + / 5 (f) • order,- • stressl., + fe(t) ■ order j • stress2j 
+ A(t) • order., • vowel.,-. 

Thus, in total, eight covariates characterize the 16 target words. 

4.3 Application results 

Our estimation yields two and three FPCs for the fRI for speakers and for the smooth 
error, respectively. No FPC is chosen for the fRI for target words. It is likely that 
the eight covariate and interaction effects describe the target words sufficiently, as 
confirmed by obtaining one FPC for the fRI for target words in the model without 
covariate effects. Most variability is explained by the curve-specific deviation which 
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also captures interactions between speakers and target words. The variance decom¬ 
position is given in Table 1. 

The left panel of Figure 2 shows the effect of covariate order (fi), which has the largest 
effect on the index trajectories. Covariate order is dummy-coded with reference cat¬ 
egory /s#sh/. Thus, the mean curves of target words with order /sh#s/are pulled 
towards the ideal reference /sh/ during the first consonant and differ slightly from the 
ideal /s/ during the second consonant compared to order /s#sh/. We conclude that 
there is an asymmetry of consonant assimilation with respect to the consonant order 
and that the assimilation is stronger for order /s#sh/. These results are consistent 
with results for English obtained by Pouplier et al. (2011). 

Moreover, we find that assimilation is stronger for target words with unstressed fi¬ 
nal target syllables (/ 2 ), for which the mean curves are pulled away from the ideal 
reference first consonant and slightly away from the ideal second consonant. This 
effect is stronger for order /s#sh/(/ 5 ). The assimilation is not affected by the stress 
of the initial target syllable for order /s^sh/(/ 3 ). For order /sh#s/(/ 6 ), however, the 
mean curves are pulled away from the ideal reference /s/ for unstressed initial target 
syllables compared to stressed initial target syllables, i.e. assimilation of the second 
consonant is increased. For both consonant orders, the vowel context mainly affects 
the transition between the two consonants (/ 4 and f-). The first consonant is closer 
to the ideal reference value in the ai compared to the ia condition, yet the second 
consonant is pulled away from its reference value. These results show that the /i/ 
vowel perturbs an adjacent consonant away from its ideal reference pattern. 

In the right panel of Figure 2, we show the effect of adding (+) and subtracting (—) a 
suitable multiple of the first FPC for speakers to the overall mean (solid line) obtained 
by setting all covariates to 0.5. The interpretation is straight forward: speakers with a 
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Table 1: Variance explained in percent of the estimated total variance. 


<7 


13.16% 7.29% 44.02% 17.11% 6.16% 8.88% 




Figure 2: Left: Effect of covariate order (red solid line) with point-wise confidence 
bands (dashed lines). Right: Mean function (solid line) and the effect of adding (+) 
and subtracting (—) a suitable multiple (2 y/bf) of the first FPC for speakers. 

negative weight for the first FPC distinguish better between the two consonants. The 
estimates for the basis weights can be used for further analysis. Further application 
results including plots for all mean effects can be found in Appendix C. 


5 Simulations 

5.1 Simulation designs 

We conduct extensive simulation studies to investigate the performance of our method. 
The data generating processes can be divided into two main groups: 1) data that mim¬ 
ics the irregularly sampled consonant assimilation data and 2) sparsely sampled data 
with a higher number of observations per grouping level but fewer observations per 
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curve. For all settings, we generate 200 data sets. 

Application-based simulation scenarios. We consider two application-based scenarios, 
one with an £RI for speakers and covariate mean effects (fRI scenario) and another 
with crossed fRIs for speakers and for target words, respectively, but no covariate 
mean effects (crossed-fRIs scenario). We generate the data based on the estimates 
of model (2.2) for our consonant assimilation data with fi(t, Xijh) corresponding to 
(4.1) and to a simple smooth intercept //(f), respectively. The data analysis yields 
two FPCs for the fRI for the speakers and three FPCs for the smooth error term. For 
the crossed-fRIs scenario, we additionally obtain one FPC for the fRI for the target 
words. The FPC weights and the measurement errors are independently drawn from 
normal distributions with zero mean and with the respective estimated variances. 
Note that in similar FPC-based approaches for simpler models/non-sparse data (e.g. 
Yao et ah, 2005; Greven et ah, 2010), simulation results were not dependent on the 
normal distribution. More details on the data generation can be found in Sections 
4.3, 5.2, and in Appendices C and D. 

Sparse simulation scenario. In order to investigate the estimation performance in the 
sparse case, we additionally generate data with crossed fRIs as in (2.2) consisting 
of observations that are sparsely sampled on [0,1]. For Bi(t) and Cj(t), we choose 
/ = J = 40 replications each with each combination observed Hij = 3 times. The 
number of observation points per curve is drawn from the discrete uniform distribution 
U{ 3,10}. For each grouping variable, we use two FPCs to generate the underlying 
process. Eigenvalues are generated as Vk = 2 A, k = 1, 2 for all three random processes. 
We choose normalized Legendre Polynomials adapted to the interval [0,1] as FPCs 
for Bi(t) and Cj(t). The first two elements of the orthonormal bases are 
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4>f (t) — 1 4\{4) = V / 3(2t — 1) 4*1 if) = V / 2sin(27rt) 

4*244) = \/5(6f 2 — 6t + 1) 4>2(4) = v^7 ( 20f 3 — 30 1 2 + 12 1 — 1) 4>2(4) = \/2cos(27r t). 


Note that the different bases need not be mutually orthogonal. The FPC weights 
and the measurement errors are independently drawn from the normal distributions 
A/”(0, Uk) and J\f(0,a 2 ), respectively. For the smooth error term, Eijh(t ), we choose a 
basis of sine and cosine functions. No covariates are included in the mean function 
pit) = sin(t) + t. We set the error variance to cr 2 = 0.05. 

For all scenarios, we center the FPC weights such that the weights of each grouping 
variable also empirically have zero mean. Moreover, we decorrelate the basis weights 
belonging to one grouping variable and assure that the empirical variance corresponds 
to the respective eigenvalue. This is done to obtain data that meets the requirements 
of our model. It allows us to separate the effect of unfavourably drawn weights and of 
the estimation performance. This adjustment gains importance for small sample sizes 
/, J, and n and also when the true eigenvalues are high. Note that in practice, we do 
not have centred and decorrelated FPC weights and thus estimates for small sample 
sizes will reflect the distribution in the sample rather than that in the population. To 
assess the impact of this procedure, we also compare our results to those of simulations 
using the original (non-centred and non-decorrelated) FPC weights, which can be 
found in Appendix D. 

We fix the number of FPCs in order to separate the effect of the truncation from the 
estimation quality. We use five marginal basis functions each for the estimation of 
the auto-covariances and eight basis functions for the estimation of the mean. We 
predict the FPC weights as EBLUPs for all scenarios and additionally compare with 
the computationally more expensive FAMM prediction (FPC-FAMM) for the fRI 
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scenario with covariates. 

We compare our FPC-based approach to a spline basis representation of the functional 
random effects (using eight basis functions) within the FAMM framework of Scheipl 
et ah (2015) (spline-FAMM). To the best of our knowledge, the work of Scheipl et al. 
(2015) implemented in the R-function pffr is the only competitor to our approach as 
all other methods are either restricted to equal, fine grids or do not allow for a crossed 
structure. Due to the high computational costs of Scheipl et al. (2015), we restrict 
our comparison to the fRI scenario, in which we can compare estimation quality and 
CBs coverage for covariate effects. 

5.2 Simulation results 

We focus our discussion on the FPC-based results for the application-based scenario 
with crossed fRIs and compare with the other settings and estimation approaches. 
We use root relative mean squared errors (rrMSE) as measures of goodness of fit 
which are of the general form yj (true-estimated) 2 /t rue 2 . rrMSE definitions for scalars, 
vectors, and bivariate functions are given in Appendix D. For the simulations of the 
fRI scenario with covariate effects, we additionally evaluate the average point-wise 
and the simultaneous coverage of the point-wise CBs. The complete results for all 
simulations are given in Appendix D. 

Simulation results for the crossed-fRIs scenario . Figure 3 shows the true and esti¬ 
mated FPCs of the two fRIs as well as of the smooth error term. As expected, the 
FPCs are estimated better the more independent levels there are for the correspond¬ 
ing grouping variable which can enter the estimation of the auto-covariance. The 
FPCs of the smooth error term (707 levels) are estimated best, followed by the FPC 
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Figure 3: True and estimated FPCs of the crossed fRIs B,(t) and C 3 (t) (top row), 
as well as of the smooth error Eijh(t ) (bottom row). Shown are the true functions 
(red solid line), the mean of the estimated functions over 200 simulation runs (black 
dashed line), the point-wise 5th and 95th percentiles of the estimated functions (blue 
dashed lines), and the estimated functions of all 200 simulation runs (grey). 


of the fRI for target words (16 levels). Most variability in the estimates is found for 
the FPCs of the fRI for speakers due to the small number of speakers (1 = 9), but 
the main features of the curves are still recovered relatively well. We obtain similar 
results for the fRI scenario. The number and complexity of the FPCs also plays an 
important role for the estimation quality, as can be seen from the results for the 
sparse scenario, where the first FPC of B t (t) (40 levels) is estimated better than the 
first FPC of Eijh(t ) (4800 levels). The latter has a more complex form, difficult to 
capture with five basis functions. 

Table 2 lists the rrMSEs averaged over 200 simulation runs for all model components. 
It shows that the mean function is reconstructed very well, which is also the case in 
the sparse scenario. The covariate effects for the fRI scenario are discussed below. 
The auto-covariances and their eigenvalues have similar low average rrMSEs for both 
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Table 2: rrMSEs averaged over 200 simulation runs for all model components by 
random process. Rows 1-3: Number of grouping levels L x and average rrMSE for 
Bi(t ), Cj(t ), and Ejjhit) and their covariance decompositions. Last row: Average 
rrMSEs for Yijh(t ), and a 2 . 


X 

L x 

K x 




v? 

V 2 

^3 




X 

(i a 2 

B 

9 

0.26 

0.15 

0.18 


0.15 

0.34 


0.18 

0.35 


0.22 


C 

16 

0.32 

0.05 



0.31 



0.12 



0.13 


E 

707 

0.06 

0.02 

0.03 

0.02 

0.04 

0.08 

0.03 

0.17 

0.19 

0.26 

0.19 


Y 












0.10 

0.02 0.09 


application-based scenarios. For the sparse scenario, the eigenvalues are estimated 
even better with average rrMSEs between 0.02 and 0.05. For the auto-covariances 
for the sparse scenario, we obtain average rrMSEs of 0.06 for each of the crossed fRIs 
and an average of 0.14 for the smooth error which is due to the complex eigenfunc¬ 
tions mentioned above. The error variance has similar low average rrMSEs for the 
two application-based scenarios. For the sparse scenario, the average rrMSE is higher 
which is due to the estimation inaccuracies in the auto-covariance of the smooth error. 
The prediction quality of the basis weights clearly depends on the estimation quality 
of the FPCs and of the eigenvalues, as well as of the error variance, as evident from 
(3.6). Also important for the prediction of the basis weights is the number of curves 
with the given weight entering the prediction. Thus, the basis weights of Cj(t ) are 
better predicted than those of E^hit). As expected, basis weights of FPCs that ex¬ 
plain more variability are estimated better. Similar results can be found for the fRI 
and for the sparse scenario. 

For all scenarios, we obtain good results for the functional random effects as well as 
for the functional response. The rrMSEs for the functional response are lowest, which 
is due to the fact that even if the FPC bases are not perfectly estimated, they can 
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still serve as a good empirical basis. Thus, the data can be reconstructed very well. 
We found considerably more outliers of the relative errors for the sparse scenario than 
for the other two scenarios, which is most probably due to an unfavourable distribu¬ 
tion of the few observation points across the curves in a few data sets. Overall, we 
can conclude that all components are estimated well and especially for the functional 
response we obtain very small rrMSEs across all simulations. 

Comparison of the different estimation results for the fRI scenario. We find that the 
functional random processes as well as the functional response are estimated equally 
well for the two options of the basis weights prediction. The functional response 
is again estimated very well with an average rrMSE of 0.09 for both EBLUP and 
FPC-FAMM estimation. The spline-FAMM results are considerably worse for the 
random processes (almost three (smooth error) and almost seven (fRI) times higher 
average rrMSEs), which results from the fact that the constraint Xi(t) = 0, X E 
{B, E}, is not fulfilled and parts are shifted between terms. The functional response is 
recovered reasonably well, but has a more than 1.5 times higher average rrMSE than 
the EBLUP and FPC-FAMM estimates. Note that due to high computation times 
(see below), we only consider 100 simulation runs for the spline-FAMM simulation. 
For the covariate effects, the FPC-FAMM estimation gives better results than the 
estimation under an independence assumption and considerably better results than 
the spline-FAMM estimation (between 2.8 and six times higher average rrMSEs). 
The average point-wise coverage of the point-wise CBs is very good for most effects 
for FPC-FAMM and the simultaneous coverage is reasonable. Both are considerably 
better than for the spline-FAMM alternative. The coverage for the latter would most 
probably improve by increasing the number of spline basis functions which is, however, 
limited by the high computation time. 
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Computation times. Our simulations show that the FPC-based approach has clear 
advantages in terms of computational complexity, despite the computational cost of 
auto-covariance estimation. We compare times for one simulation run of the fRI 
scenario for each estimation option obtained under the same conditions (without 
parallelization in function bam that would speed up the estimation). The study was 
run on a 64 Bit Linux platform with 660 Gb of RAM memory. The FPC-based 
approach with the basis weights predicted as EBLUPs took five hours and predicting 
the basis weights using FPC-FAMM took slightly less than ten hours longer. The 
spline-FAMM took by far the longest with a duration of ten days which is due to the 
two extra smoothing parameters each for the fRI and the smooth error which have to 
be estimated. Moreover, using FPCs reduces the number of necessary basis functions. 


6 Discussion and outlook 

We propose an FPC-based estimation approach for functional linear mixed models 
that is particularly suited to irregularly or sparsely sampled observations. To pool 
information, we smooth both the mean and auto-covariance functions. We propose 
and compare two options for the prediction of the FPC weights and obtain conditional 
point-wise confidence bands for the functional covariate effects. Our simulations show 
that our method reliably recovers the features of interest. The parsimonious repre¬ 
sentation of the functional random effects in bases of eigenfunctions outperforms the 
spline-based alternative of Scheipl et al. (2015) with which we compare, both in terms 
of error rates and coverage as well as in terms of computation time. To the best of 
our knowledge, there is no other competitor to our approach as all other methods are 
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either restricted to regular grid data or simpler correlation structures. In our applica¬ 
tion to speech production data, we show that our method allows to draw conclusions 
about the asymmetry of consonant assimilation to an extent which is not achievable 
using conventional methods with data reduction. 

Building on existing methods for our estimation approach allows us to take advantage 
of robust, flexible algorithms with a high functionality. The computational efficiency, 
however, could potentially be improved by exploiting the special structure of our 
model. In future work, we plan to improve the estimation of the auto-covariances in 
order to better account for their symmetry and positive semi-definiteness and for the 
fact that the cross products in (3.3) are not homoscedastic. Moreover, it would be 
interesting to compare the different options for iterative estimation in detail. 

The construction of point-wise and simultaneous confidence bands that account for 
the variability of the estimated FPC decomposition is beyond the scope of this work, 
but would be of interest. For uncorrelated functions, Goldsmith et al. (2013) propose 
bootstrap-based corrected confidence bands for densely and sparsely sampled func¬ 
tional data. However, it remains an open question how to extend their non-parametric 
bootstrap to our correlated curves, and computational cost is another issue. 
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A Derivations 


Derivation of the variance decomposition of model (2.2). Using iterated expectations 
allows us to decompose the variance of the response as follows 


Jvai[Y ijh (t)]dt = J^Yar[B i (t)]dt +J^av[C j (t)]dt +J^Yav[E ijh (t)]dt 


'T 

+ J £ ijh (t)dt 

OO « 


< i>k( t )<i>k( t ) dt +' 52 i/ k / $(t)<t>k(t)dt 

k =1 Jr ■ *~i 


k =1 


=1 


=1 


OO /» p 

+ / <f>k( t )<f>k( t ) dt + / (j2dt 

fc=l _, Jr 


= 1 

OO OO 


= E-f+Eif+E-f+^in 


k =1 /c=l /c=l 


A.l Marginal bases for the estimation of the mean 

^f p g ,p = l,...,P, is an inflated vector of length £> containing the values of covariate 
x p = i i p ,.... XuHjjp} 7 ■ For each curve, the covariate values are replicated for 
each observation point yielding 

%lllp 

^lllp 

\ifP — : 

9 

XlJHjjp 


XlJHup 

For the functional intercept fo(t), = 1®. 


p = 0,... ,P, is aD x K p matrix comprising the evaluations of the basis functions 
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on the original observations points of the form 


V’i(bm) 

V’Kp(bui) 

V’llbiiTm) 

bxp(bllTm) 

^i(bjffjji) 

Vkp^IJHui) 

V’i (tiJHuTumj) •• 

■ Vkp^IJHuTuh^) 


A. 2 Matrices in the prediction of the basis weights as EBLUPs 


is a block-diagonal matrix of dimension T> x IN 1 


= 


4>f (bm) 


^B(biii) 


h (tl JH 1 jT 1 jh 1 j) 


(bin) 


^wB(bin) 


01 (0 IJHijTijhjj ) ( t ) N B (t'IJHijTijH IJ ) 
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<fr c consists of / block-diagonal matrices - one for each speaker - and is of dimension 
D x JN C . 


$ c = 


4 >\ (tin) 


4 >? (bin) 




(tllffuTiiHn ) $N C (bmuTiiH!! ) 


(tmi) 


(>l(t I1HllTllHii ) ■■■ 


4>1 (tun) 


1 


<t>\ (tun) 


djv c (tun) 


djyc (tun) 


L 4>l(tlJH IJTlJHij ) ••• 4> C N c(tlJH IJTlJHij ) 

Obviously, the role of speakers and target words is exchangeable. For present pur¬ 
poses, however, we assume that all vectors and matrices are first ordered by speakers 
and within each speaker ordered by target words. is a block-diagonal matrix of 
dimension £> x nN E with blocks 


df(bm) ••• df E (hm) 


1 

C~t- 

tq 

tq 

c-+* 

•-I 

tq 

_1 

_ ^ffilllTm) ••• 4 *% e (tlllTm ) 

} # * # } 

_ df (tlJHuT IJHlJ ) ••• ^EdlJHuTuujj) _ 


Note that for irregularly spaced functional data, the FPCs evaluated on the original 
observation points are curve-specific. 

The estimated covariance matrix of the £, G, is of the form 


G = 


Cov(0 

C3v(£ c ) 

C3v(£ e ) 
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where each covariance matrix Cov(£ x ),X G {B,C,E} is a diagonal matrix with 
elements z>f,..., v* x ,..., z>f,..., v* x . 

A. 3 Marginal bases for FPC-FAMM 

d/^ is a D x / incidence matrix of which the entries in the cJtli column are one wherever 
the row belongs to the c/th speaker and zero otherwise, d/^ and d/^ are analogously 
D x J and D x n matrices with entries in the c/tli column equal to one wherever the 
row belongs to the citli target word or to the c/th curve, respectively. For lack of space, 
we exemplarily show in the following the part of the matrices d/^, X G {B,C,E}, 
that belongs to the ith speaker (denoted by d/^ ? j, assuming that the data are ordered 
by speakers, within each speaker ordered by target words, within each target word 
ordered by repetition, then by curves, and finally by observation points. 


r 

l 

_ 




_ 




_ 

0 ••• 

1 

... 0 


1 

0 ••• 

0 


1 

0 ••• 

0 

0 ••• 

1 

... 0 


1 

0 ••• 

0 


1 

0 ••• 

0 

0 ••• 

1 

... 0 


1 

0 ••• 

0 


0 

0 ••• 

0 

0 ••• 

1 

... 0 

\D Ci — 

9 

1 

0 ••• 

0 

II 

* 

0 

0 ••• 

0 

0 ••• 

1 

... 0 


0 

0 ••• 

1 


0 

0 ••• 

0 

0 ••• 

1 

... 0 


0 

0 ••• 

1 


0 

0 ••• 

0 

0 ••• 

1 

... 0 


0 

0 ••• 

1 


0 

0 ••• 

1 

0 ••• 

1 

... 0 


0 

0 ••• 

1 


0 

0 ••• 

1 


d/^ then contains the stacked partial matrices for all speakers. 
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is a © x N x matrix containing the evaluations of the respective eigenfunctions 
on the original observation points of the form 


= 


<f>i (bin) 


<I>1 (bllTm) 


4>l (tijHu l) 


^x(tllll) 


^x(bllTiu) 




(tuHjjT, JH ) ••• (tuHuTj JH ) 


, Xe{B,C,E}. 


B Supplementary details on the estimation and imple¬ 
mentation 

B.l Implementation of the auto-covariance estimation 

Model (3.3) does not contain an intercept. The implementation for the fRI design, 
however, requires that an intercept is included and added to the auto-covariance of 
the fRI due to the centring constraint in function bam. This is not the case for the 
crossed-fRI design, where each smooth is varied by an indicator variable and thus no 
constraint is applied. For more details see the description of the function gam and of 
gam.models. 
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B.2 Rescaling of the eigenvectors and eigenvalues 


>) 


T 


Denote the estimated eigenvectors by <j>£ = 




, and the estimated 


eigenvalues by , k G N, X £ { B , C, E\. In order to ensure that the approximated 
functions are orthonormal with respect to the scalar product (f,g) = f f(t)g(t)dt, 
we rescale the eigenvectors by 



with a denoting the constant interval width of the equidistant grid (fi,..., tp}. As 
a consequence, the eigenvalues to the rescaled eigenvectors need to be adjusted as 



B.3 Truncation of the FPCs 

Due to the additive structure in the variance decomposition, we can choose the trun¬ 
cation lags N B , N c , and N E in the following way: 

1. specify the proportion of explained variance L, e.g. L = 0.95 as used in the 
application in Section 4. 

2. select the FPCs corresponding to their eigenvalues in decreasing order, until 



(B.l) 


Note that the three random processes are treated equally, i.e. in each step, the FPC 
with the highest corresponding eigenvalue is chosen regardless of the associated pro¬ 
cess. In practice, all infinite sums in (B.l) are approximated by the finite sums of all 
obtained eigenvalues. The eigenvalues and the error variance are replaced by their 


estimates. 
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B.4 Fixing the smoothing parameter in FPC-FAMM 

Technically, in order that a fixed smoothing parameter A* (here A* = 1) is used in 
the prediction of the FPC-based fRIs, we need to specify the smoothing parameter 
in function pffr as 

A fi x = A* • S.scale ■ d 2 , (B.2) 

where S.scale is a scaling factor used in the set-up of the design matrices in order to 
numerically stabilize the prediction. It can be obtained by setting f it=FALSE in the 
call of function pffr. d 2 is an estimate of the error variance. We can use the estimate 
obtained in step 2 of our estimation procedure. Note that equation (B.2) makes clear 
that the point-wise confidence bands for the mean function are not only conditional 
on the estimated FPCs and the truncation level, but also on the estimated error 
variance. 

B.5 Iterative estimation 

If desired, the estimation accuracy can be improved by applying the estimation steps 
iteratively. At least two possibilities exist: We can either perform steps 1-4 and then 
re-start with step 1 with the adjusted observations Y*- h {t) := Yijhif) ~ Bi(t) — Cj(t ) — 
Eijh(t ) until a pre-defined criterion is reached. Alternatively, we can replace the mean 
with that obtained in the FAMM framework in step 4 and restart with step 2. 
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C Supplementary application details and results 

C.l Pre-processing 

The index calculation is based on the calculation of the power spectrum over a time 
window of approximately 20 ms, shifted in 5 ms steps over the time interval of the 
consonants. In order to be able to compare the index curves for the two consonant 
orders (/s#sh/, /sh#s/) directly, the index curves of /sh^s/were mirrored along the 
time axis (mapping +1 to -1 and vice versa) such that for both orders the index 
dynamic ranges from +1 for the first consonant to -1 for the second consonant rather 
than from +1 for /s/ to -1 for /sh/. To achieve this, we first estimated smooth mean 
curves of available reference curves of orders /sh#sh/ and /s#s/ per speaker and for 
each combination of the covariates vowel, stressl, stress2 using penalized splines and 
evaluated them on the measurement points. We then mirrored the index curves of 
order /sh#s/at the speaker-condition specific mean values, averaging over the mean 
curves for /sh^s/and /s#sh/. 

C.2 Supplementary application results 

In the following, we show additional application results, including the eigenvalues, 
the effects of the covariates and interactions, the second FPC for speakers, and the 
variance decomposition of the model without covariates on which we base our crossed- 
fRI simulation setting. 

Table 3 gives the estimated eigenvalues for the model with covariates included. Figure 
4 shows the estimated effects and point-wise confidence bands of the covariates stressl 
(0:Strong, l:Weak), stress2 (0:Strong, l:Weak), and vowel (0: ia, l:ai) as well as of 
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Figure 4: Covariate mean effects (red solid lines) with conditional point-wise con¬ 
fidence bands (dashed lines). Upper: Reference mean fo(t). Middle (from left to 
right): covariate effects of covariates stressl, stress2, vowel. Lower (from left to 
right): interaction effects of order and stressl, order and stress2, order and vowel. 

the interactions between covariate order and the other three covariates. 

The second FPC for speakers is depicted in Figure 5. The interpretation is as follows: 
The index curves of speakers with positive FPC weights for the second FPC are pulled 
towards the ideal reference value of the second consonant, whereas the index curves 
of speakers with negative FPC weights are pulled towards the ideal reference value of 
the first consonant. 

Table 4 gives the variance decomposition of the model with no covariate included. 
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Table 3: Estimated eigenvalues for the model with covariates. 




<J 


5.84 TO- 3 3.23T0“ d 19.53T0“ d 7.59T0“ d 2.73T0" 3 3.94 TO 


1-3 


1-3 


-3 


1-3 


1-3 



Figure 5: Second FPC for speakers. Shown is the mean function (solid line) and the 
effect of adding (+) and subtracting (-) a suitable multiple (2y / hjf) of the first FPC 
for speakers. 


Table 4: First row: Estimated eigenvalues of K B and K E and estimated error vari¬ 
ance for the crossed model without covariates. Second row: Variance explained in 
percent of the estimated total variance. 


'C 


<7 


5.86-IQ- 3 2.71-IQ- 3 8.89 • 10“ 3 19.1 • 10“ 3 7.53 • 10“ 3 2.66 • 10“ 3 5.62 • 10“ 3 


10.83% 


5.00% 16.44% 35.23% 13.93% 


4.92% 10.39% 
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D Supplementary simulation details and results 


D.l Measures of goodness of fit 

We use the root relative mean squared error 



(D.l) 


as a measure of goodness of fit for vector-valued estimates 0 of 6 = (#i,... ,$l) T 
(FPC weights Ck = (£&, • • •, , k = 1,, N x , X G {B, C, E}) and for scalar 


estimates (eigenvalues is*, k = 1 ,..., N x , X e {D, C, E} and the error variance 
a 2 ) as a special case with L = 1 . For the FPC weights, (D.l) is approximately 



For all functions 9{t) (mean function with and without covariates fo(t ),..., / 7 (f), 
eigenfunctions (f>k(t), k = 1 ,..., X e {D, 6 ', F 1 }), we approximate the integrals 
by sums and obtain 




(D.2) 




Note that for the eigenfunctions, the denominator simplifies to approximately one 
as f 4>k(t) 2 df = 1 . As the eigenfunctions are only unique up to sign, we mirror 
the estimated eigenfunctions around the x-axis, also compute the rrMSE for the 
mirrored estimates, and choose the smaller rrMSE. For the fRIs and for the response 
function, we additionally average over the respective levels. As the fRIs are centred, 
the denominator simplifies to the average variance. The rrMSE for the bivariate 
functions (auto-covariances) are defined analogously. 
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D.2 Results for simulations with centered and decorrelated 
basis weights 


Simulation results for the crossed-fRIs scenario . In the following, additional simula¬ 
tion results for the application-based crossed-fRIs scenario with centred and decorre¬ 
lated basis weights are shown. Figure 6 shows the true and estimated mean functions. 
In Figure 7, the boxplots of the estimated eigenvalues for each Cj(t), and 

are depicted. In Figure 8, we show the boxplot of the estimated error variances. 



Figure 6: True and estimated mean functions. Shown are the true function (red), 
the mean of the estimated functions over 200 simulation runs (black dashed line), the 
point-wise 5th and 95th percentiles of the estimated functions (blue dashed lines), 
and the estimated functions of all 200 simulation runs (grey). 




48 Jona Cederbaum et a 1. 



Figure 7: Boxplots of the estimated eigenvalues of the auto-covariances of the crossed 
fRIs (top row), as well as the eigenvalues of the auto-covariance of the smooth error 
(bottom row) for all 200 simulations runs. 



Figure 8: Boxplots of the estimated error variances a 2 for all 200 simulation runs. 
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Simulation results for the fRI scenario . In the following, we show additional simu¬ 
lation results for the application-based fRI scenario with centred and decorrelated 
basis weights. Figure 9, Figure 10, and Figure 11 show the true and estimated co¬ 
variate and interaction effects based on the independence assumption (Figure 9), on 
FPC-FAMM estimation (Figure 10), and on the spline-FAMM alternative (Figure 
11), respectively. We evaluate the performance of the point-wise CBs obtained by 
FPC-FAMM and spline-FAMM, by looking at the point-wise coverage shown in Fig¬ 
ures 12 and 13. We additionally compare the simultaneous coverage of the point-wise 
CBs in terms of percentage of completely covered curves in Table 5. 

Figure 14 depicts the true and estimated FPCs of the fRI and of the smooth error 
and Figure 15 shows the boxplot of the estimated eigenvalues. The boxplot of the 
estimated error variances is shown in Figure 16. Table 6 lists the average rrMSEs for 
all model components except for the covariate effects and Table 7 lists the rrMSEs 


for the covariate effects. 
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Figure 9: True and estimated covariate and interaction effects estimated using the 
independence assumption. Shown are the true function (red), the mean of the esti¬ 
mated functions over 200 simulation runs (black dashed line), the point-wise 5th and 
95th percentiles of the estimated functions (blue dashed lines), and the estimated 
functions of all 200 simulation runs (grey). 
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Figure 10: True and estimated covariate and interaction effects estimated using FPC- 
FAMM. Shown are the true function (red), the mean of the estimated functions over 
200 simulation runs (black dashed line), the point-wise 5th and 95th percentiles of 
the estimated functions (blue dashed lines), and the estimated functions of all 200 
simulation runs (grey). 
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Figure 11: True and estimated covariate and interaction effects estimated using spline- 
FAMM. Shown are the true function (red), the mean of the estimated functions over 
100 simulation runs (black dashed line), the point-wise 5th and 95th percentiles of 
the estimated functions (blue dashed lines), and the estimated functions of all 100 
simulation runs (grey). 
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Figure 12: Average point-wise coverage of the point-wise CBs obtained by FPC- 
FAMM for all covariate and interaction effects. For each effect, the point-wise cover¬ 
age averaged over 200 simulation runs (black line) is shown. The red line indicates 


the nominal value of 0.95. 
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Figure 13: Average point-wise coverage of the point-wise CBs obtained by spline- 
FAMM for all covariate and interaction effects. For each effect, the point-wise cover¬ 
age averaged over 100 simulation runs (black line) is shown. The red line indicates 


the nominal value of 0.95. 
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Table 5: Simultaneous coverage of the point-wise CBs for FPC-FAMM and for spline- 

FAMM. Shown is the proportion of completely covered curves for all covariate and 

interaction effects. For FPC-FAMM, the coverage refers to 200 simulation runs, 

whereas for splines-FAMM, 100 simulation runs are taken into account. 

fo(t) fi(t) h(t) f 3 (t) f 4 (t) f 5 (t ) f 6 (t) f 7 (t ) 

n(t, x ijh ) FPC-FAMM 43.5% 71.5% 70.5% 64.5% 76.0% 70.0% 66.0% 78.5% 


/i(t, Xijh ) spline—FAMM 


1.0% 5.0% 1.0% 1.0% 2.0% 2.0% 2.0% 5.0% 






Figure 14: True and estimated FPCs of the fRI (top row) and of the smooth error 
(bottom row). Shown are the true functions (red), the mean of the estimated functions 
over 200 simulation runs (black dashed line), the point-wise 5th and 95th percentiles 
of the estimated functions (blue dashed lines), and the estimated functions of all 200 
simulation runs (grey). 
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Figure 15: Boxplots of the estimated eigenvalues of the auto-covariances of the fRI 
(top row), and of the smooth error (bottom row) for all 200 simulations runs. 



Figure 16: Boxplot of the estimated error variances a 2 for all 200 simulation runs. 

Table 6: rrMSEs averaged over 200 simulation runs for all model components by 
random process. Rows 1-3: Number of grouping levels L x and average rrMSEs for 
the fRI and the smooth error. Last row: Average rrMSEs for the functional response, 
the mean, and the error variance. 
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0.02 

0.02 
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Table 7: Average rrMSEs for the estimated mean and covariate effects. Rows 1- 
2: rrMSEs for the estimation based on the independence assumption and for the 
estimation using FPC-FAMM averaged over 200 simulation runs. Last row: rrMSEs 
for the spline-based alternative averaged over 100 simulation runs. 



/o(t) h(t) f 2 (t) hit) U(t) f 5 (t ) f 6 (t) f 7 (t ) 

/i(t, Xijf 

0.06 

0.13 

0.22 

1.43 

0.30 

0.58 

0.39 

0.60 

M t , ^zjTi) FPC-FAMM 

0.06 

0.12 

0.20 

1.34 

0.25 

0.53 

0.39 

0.47 

/i(t, spline-FAMM 

0.36 

0.38 

0.58 

3.83 

0.84 

1.54 

1.10 

1.85 


The average rrlSEs of the covariate and interaction effects lie between 0.06 {fo{t)) 
for both estimation options and 1.43 or 1.34 (/ 3 (t)) for the estimation using the 
independence assumption or FPC-FAMM, respectively. Note that the high value for 
f 3 (t) is the result of the fact that the true covariate effect is very close to zero along the 
whole time interval, and to avoid dividing by values near zero it is more meaningful 
to look at the root mean squared error (rMSE) instead which is similar to the rMSEs 


of other covariates. 
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Simulation results for the sparse scenario. In the following, additional results for the 
sparse scenario with centred and decorrelated basis weights are shown. Figure 17 
shows the true and estimated mean functions and Figure 19 depicts the boxplot of 
the estimated eigenvalues for the two fRIs and for the smooth error. In Figure 20, we 
show the boxplot of the estimated error variances. In Table 8, the average relative 
errors for all model components are given. 



Figure 17: True and estimated mean function pit, a;^). Shown are the true function 
(red), the mean of the estimated functions over 200 simulation runs (black dashed 
line), the point-wise 5th and 95th percentiles of the estimated functions (blue dashed 
lines), and the estimated functions of all 200 simulation runs (grey). 


Table 8: rrMSEs averaged over 200 simulation runs for all model components by 
random process. Rows 1-3: Number of grouping levels L x and average relative errors 
for the functional random effects and their covariance decompositions. Last row: 
average rrMSE of the functional response, the mean, and the error variance. 


X 

L x 

K x 4> x u x v x i x i x X 

/L (J 2 

B 

40 

0.06 0.05 0.07 0.02 0.04 0.04 0.11 0.06 


C 

40 

0.06 0.07 0.11 0.03 0.05 0.23 0.25 0.21 


E 

4800 

0.14 0.11 0.07 0.02 0.05 0.30 0.19 0.29 


Y 


0.09 

0.03 1.81 
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Figure 18: True and estimated FPCs of the crossed fRIs Bj(t) (top row) and Cj(t) 
(middle row), as well as the FPCs of the smooth error Eijh(t ) (bottom row). Shown 
are the true functions (red), the mean of the estimated functions over 200 simulation 
runs (black dashed line), the point-wise 5th and 95th percentiles of the estimated 
functions (blue dashed lines), and the estimated functions of all 200 simulation runs 
(grey). 
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Figure 19: Boxplots of the estimated eigenvalues of the auto-covariances of the crossed 
fRIs Bjt) (top row), Cj{t) (middle row), as well as the eigenvalues of the auto¬ 
covariance of the smooth error Eijh(t ) (bottom row) for all 200 simulations runs. 



Figure 20: Boxplots of the estimated error variances a 2 for all 200 simulation runs. 
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D.3 Results for simulations with original basis weights 

Before we show the results for the simulation with the original (non-centred and non- 
decorrelated) basis weights in detail, we give a short summary and compare with the 
results for the simulations with centred and decorrelated basis weights. 

We observe that with the original basis weights, the estimated FPCs are often per¬ 
muted within one grouping variable, e.g. the first and second FPC of the speakers are 
interchanged or are linear combinations of them, due to the correlation in the basis 
weights. As expected, this effect increases the smaller the corresponding number of 
independent levels, as the empirical FPC weights then have an empirical distribution 
far from the theoretical one. Moreover, we obtain higher average rrMSEs for the 
eigenvalues when using the original basis weights. Correspondingly, we also obtain 
worse results for the auto-covariances. For all three simulation settings, the rrMSEs 
for the basis weights tend to be higher with the original basis weights. Using the 
original basis weights results in higher rrMSEs for the functional random effects, es¬ 
pecially for functional random effects with a small number of grouping levels. The 
estimation of the mean function for non-centred basis weights is much worse due to a 
shift of the mean by the respective FPC multiplied by the empirical mean of the basis 
weights. Correspondingly, also the coverage of the point-wise CBs decreases for most 
points. Yet, the functional response for the original basis weights is again estimated 
very well as shifts between mean and non-centred random effects cancel out. For the 
covariate effects and the coverage of the CBs, we hardly observe any differences to the 
simulations with centred and decorrelated basis weights. Also for the error variance, 
we do not observe a considerable change in the rrMSEs. 


Simulation results for crossed-fRIs scenario . In the following, we show the results 
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for the simulations of the application-based crossed-fRIs scenario with non-centred 
and non-decorrelated basis weights. Figure 21 shows the true and estimated mean 
functions and Figure 23 depicts the boxplot of the estimated eigenvalues for the two 
fRIs and for the smooth error. In Figure 24, we show the boxplot of the estimated 
error variances. In Table 9, the average rrMSEs for all model components are given. 



Figure 21: True and estimated mean function p(t,Xijh)- Shown is the true function 
(red), the mean of the estimated functions over 200 simulation runs (black dashed 
line), the point-wise 5th and 95th percentiles of the estimated functions (blue dashed 
lines), and the estimated functions of all 200 simulation runs (grey). 


Table 9: rrMSEs averaged over 200 simulation runs for all model components by 
random process. Rows 1-3: Number of grouping levels L x and average rrMSE for the 
functional random effects. Last row: Average rrMSE for the functional response, the 
mean, and the error variance. 


X 

L x 

K x <j> x v x u x v* & g $ X 

/i a 2 

B 

9 

0.55 0.37 0.46 0.38 0.53 0.41 0.70 0.38 


C 

16 

0.32 0.05 0.31 0.24 0.25 


E 

707 

0.09 0.04 0.05 0.04 0.06 0.09 0.05 0.19 0.21 0.27 0.20 


Y 


0.10 

0.10 0.10 
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Figure 22: True and estimated FPCs of the crossed fRIs B,(t) and Cj(t) (top row), 
as well as the FPCs of the smooth error Eijh(t ) (bottom row). Shown are the true 
functions (red), the mean of the estimated functions over 200 simulation runs (black 
dashed line), the point-wise 5th and 95th percentiles of the estimated functions (blue 
dashed lines), and the estimated functions of all 200 simulation runs (grey). 
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Figure 23: Boxplots of the estimated eigenvalues of the auto-covariances of the crossed 
fRIs Bj(t) and Cj(t ) (top row), as well as the eigenvalues of the auto-covariance of 
the smooth error Eijh(t ) (bottom row) for all 200 simulations runs. 



Figure 24: Boxplots of the estimated error variances a 2 for all 200 simulation runs. 
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Simulation results for the fRI scenario . In the following, we show additional simula¬ 
tion results for the simulations of the application-based fRI scenario with non-centred 
and non-decorrelated basis weights. Figure 25 and Figure 26 show the true and esti¬ 
mated covariate and interaction effects estimated based on the independence assump¬ 
tion, and on FPC-FAMM, respectively. We evaluate the performance of the point- 
wise CBs by looking at the point-wise coverage shown in Figure 27. We additionally 
look at the simultaneous coverage of the point-wise CBs in terms of percentage of 
completely covered curves in Table 10. 



Figure 25: True and estimated covariate and interaction effects estimated based on 
the independence assumption. Shown are the true function (red), the mean of the 
estimated functions over 200 simulation runs (black dashed line), the point-wise 5th 
and 95th percentiles of the estimated functions (blue dashed lines), and the estimated 
functions of all 200 simulation runs (grey). 














sparse FLMM 65 



Figure 26: True and estimated covariate and interaction effects using FPC-FAMM. 
Shown are the true function (red), the mean of the estimated functions over 200 
simulation runs (black dashed line), the point-wise 5th and 95th percentiles of the 
estimated functions (blue dashed lines), and the estimated functions of all 200 simu¬ 
lation runs (grey). 


Table 10: Simultaneous coverage of the point-wise CBs obtained by FPC-FAMM. 
Shown is the proportion of completely covered curves for all covariate and interaction 
effects. The coverage refers to 200 simulation runs. 



fo(t) fi(t) f 2 (t ) f 3 (t ) / 4 (i) / 5 (t) /e(t) hit) 

Mb *pfe)FPC-FAMM 

11% 71% 71.5% 65.5% 74.5% 69% 67% 78% 


Table 11: rrMSEs averaged over 200 simulation runs for all model components by 
random process. Rows 1-3: Number of grouping levels L x and average rrMSEs for 
the fRI and the smooth error. Last row: Average rrMSEs for the functional response, 
the mean, and the error variance. 


X 

L x 

K X (f) X 4> X (f) X V X U X U X £ X £ X £ X X Xppc-FAMM 

<7 2 

B 

9 

0.51 0.49 0.51 0.31 0.44 0.50 0.74 0.35 0.35 


E 

707 

0.07 0.04 0.05 0.04 0.05 0.05 0.05 0.16 0.20 0.25 0.18 0.17 


Y 


0.09 0.09 

0.11 
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Figure 27: Average point-wise coverage of the point-wise CBs obtained by FPC- 
FAMM for all covariate and interaction effects. For each effect, the point-wise cover¬ 
age averaged over 200 simulation runs (black line) is shown. The red line indicates 
the nominal value of 0.95. 


Table 12: Average rrMSEs for the estimated mean and covariate effects for the esti¬ 
mation using the independence assumption (first row) and using FPC-FAMM (second 
row) averaged over 200 simulation runs. 



/o(<) fi(t) / 2 (f) h{t) U(t) f 5 (t) f 6 (t) f 7 (t) 

h) 

0.13 

0.13 

0.22 

1.43 

0.30 

0.58 

0.39 

0.59 

nitxijh) FPC-FAMM 

0.12 

0.12 

0.21 

1.34 

0.25 

0.53 

0.38 

0.47 
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Figure 28: True and estimated FPCs of the fRI B,(t) (top row) and of the smooth 
error (bottom row). Shown are the true functions (red), the mean of the estimated 
functions over 200 simulation runs (black dashed line), the point-wise 5th and 95th 
percentiles of the estimated functions (blue dashed lines), and the estimated functions 
of all 200 simulation runs (grey). 



Figure 29: Boxplots of the estimated eigenvalues of the auto-covariances of the fRI 
Bi(t) (top row), and of the smooth error (bottom row) for all 200 simulations runs. 
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Figure 30: Boxplots of the estimated error variances a 2 for all 200 simulation runs. 
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Simulation results for the sparse scenario. In the following, additional results for the 
simulations of the sparse scenario with non-centred and non-decorrelated basis weights 
are shown. Figure 31 shows the true and estimated mean functions and Figure 33 
depicts the boxplot of the estimated eigenvalues for the two fRIs and for the smooth 
error. In Figure 34, we show the boxplot of the estimated error variances. In Table 
13, the average rrMSEs for all model components are given. 



Figure 31: True and estimated mean function pit, a;^). Shown are the true function 
(red), the mean of the estimated functions over 200 simulation runs (black dashed 
line), the point-wise 5th and 95th percentiles of the estimated functions (blue dashed 
lines), and the estimated functions of all 200 simulation runs (grey). 


Table 13: rrMSEs averaged over 200 simulation runs for all model components by 
random process. Rows 1-3: Number of grouping levels L x and average rrMSEs for 
the random processes. Last row: Average rrMSEs for the functional response, the 
mean, and the error variance. 


X 

L x 

K x u x v x & & X 

fl (J 2 

B 

40 

0.25 0.22 0.22 0.17 0.18 0.21 0.36 0.15 


C 

40 

0.24 0.24 0.26 0.17 0.20 0.36 0.48 0.28 


E 

4800 

0.14 0.11 0.07 0.03 0.05 0.30 0.20 0.30 


Y 


0.09 

0.32 1.76 
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Figure 32: True and estimated FPCs of the crossed fRIs Bj(t) (top row) and Cj(t) 
(middle row), as well as the FPCs of the smooth error Eijh(t ) (bottom row). Shown 
are the true functions (red), the mean of the estimated functions over 200 simulation 
runs (black dashed line), the point-wise 5th and 95th percentiles of the estimated 
functions (blue dashed lines), and the estimated functions of all 200 simulation runs 
(grey). 
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Figure 33: Boxplots of the estimated eigenvalues of the auto-covariances of the crossed 
fRIs Bjt) (top row), Cj{t) (middle row), as well as the eigenvalues of the auto¬ 
covariance of the smooth error Eijh(t ) (bottom row) for all 200 simulations runs. 



Figure 34: Boxplots of the estimated error variances a 2 for all 200 simulation runs. 


















































